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abstract 

The a ncestral selection graph in population genetics was introduced by 



Krone and Neuhauser 



(jl997l ) as an analogue of the coalescent genealogy of a sample of genes from a neutrally 
evolving population. The number of particles in this graph, followed backwards in time, 
is a birth and death process with quadratic death and linear birth rates. In this paper 
an explicit form of the probability distribution of the number of particles is obtained by 
using the dens ity of the allele frequency in the corresponding diffusion model obtained by 



Kimural (|l955l ) . It is shown that the process of fixation of the allele in the diffusion model 
corresponds to convergence of the ancestral process to its stationary measure. The time 
to fixation of the allele conditional on fixation is studied in terms of the ancestral process. 

Keywords: Ancestral Process, Diffusion Model, Ancestral Selection Graph, Coalescent 
Model, Birth and Death Process 



1. Introduction 



In population genetics, the genealogy of a sample of genes plays an important role in a 
probabilistic description of the sample. Consider a discrete-time Wright-Fisher model of 
a population consisting of 2N neutral genes. If we measure time in units of 2N genera- 
tions and let A^ tend to infinity, then the Wright-Fisher process converges to a diffusion 
process. The genealogy of a sample of n ge nes, followed backwards in time, is described 



by the coalescent process (JKingmanl . 



1982bl ). The convergence is robust under a number 
of different models (e.g. Moran model). Let {an{t);t > 0} be the number of ancestors of 
a sample of n genes at time t in the past. Then the size process {a„(t); i > 0} is a death 
process with death rate i{i — l)/2 when the size is i. The size pr ocess wil l be referred to 



as the ancestral process. The distribution of a„(t) is known (see 



Griffith^ (ll979|) 



Tavare 



(119M)) and 

(1.1) PK(t) = i] 



E 



{-iY-^{2k-i){i)u-Mk.,.. 



l,2,...,n. 



k=i 



i\{k-i)\{n)k 

where p\{t) := exp{—k{k — l)t/2}, [n]^ = n{n — l) ■ ■ ■ (n — k + l) and (n)^ = n{n + l) ■ ■ ■ {n-\- 
k — 1). The total variation norm between an{t) and an{oo) = Si has a simple form 



(1.2) 



\an{t),Si\ 



l-PK(t) = l]. 



There is a first time W^^ such that an{W^i) = 1, which is the time to the most recent 
common ancestor. The density of W^ i follows 

(_l)/c-i(2fc-l)[n], „ 



:i.3) 



p«i < t] = PK(t) = 1] = j; 



A bound for P[a„(t) = 1] is known (see 
(1.4) 



-k-i- 



Kingmanl (|l982bl ). 



(")a 



-pUt). 



Tavard (Il984l ^1 and 



rt 1 

e-t < 1 _ p[a^(t) = 1] < 3^— e-*, 

n + 1 



n = 2,3,... 



Krone and Neuhauseij (jl997l ) is an analogue 



The ancestral selection graph introduced by[ 
of the coalescent genealogy. Assume that a pair of allelic types Ai and A2 are segregating 
in a population, and the selective advantage of a type Ai gene over a type A2 gene is 



s (> 0). Let A^ — > oo while c = Ns is held constant. The elements are referred to as 
particles. Let bn{t) be the number of edges, or ancestral particles, in a cross section of the 
ancestral selection graph for a sample of n genes at time t in the past. In the ancestral 
selection graph, coalescence occurs at rate ai = i{i — l)/2 and branching occurs at rate 
Pi = 2ci when the size is i. Then the ancestral process {6„(t);t > 0} is a birth and death 
process with rates Pi and Qj. A particle is called real if it is a part of the real genealogy 
of the sample, otherwise the particle is called virtual. If two particles reach a coalescing 
point, the resulting particle is real if and only if at least one of the two particles is real, 
otherwise the resulting particle is virtual. If a real particle reaches a branching point, it 
splits into a real particle and into a virtual particle. If a virtual particle reaches a branching 
point, it splits into two virtual particles. If a type A2 particle reaches a branching point, it 
splits into two type A2 particles. If a type Ai particle reaches a branching point, it splits 
into two particles, where at least one of the two particles is type Ai. Because the death 
rates are quadratic while the the birth rates are on ly linear, there is a first time W^ ^ such 



that bniWfi^) 



1. 



Krone and Neuhaused (119971 ) consider stopping the process at this 



time, since the genetic composition of the sample is determined by then. They called the 
ancestral particle at the time the ultimate ancestor. In the case of no mutation, the real 
genealogy embedded iii an aii cestral selection graph is the same as in the neutral process 



(JKrone and Neuhauseii (|l997l ): Theorem 3.12). Thus the ancestral process of the real 
particles can be described by the neutral process {a„(t); t > 0}. In this article, we discuss 
properties of the ancestral process {bn (t)\t > 0} withou t mutation which is not stopped 



upon reaching the ultimate ancestor. 



FearnheadI (|2002l ) has studied a process which is 



not stopped upon reaching the ultimate ancestor. He identifies the stationary distribution 
of the process and uses the distribution to characterize the substitution process to the 
common ancestor. 



Kimural (|l955l ) studied the density of the allele frequency by analyzing the diffusion 



process to which the Wright-Fisher model with directional selection converges. Let Xp{t) 



5 

be the frequency of the ahele Ai at time t forward in a population in which the initial 
frequency is p. Then the Kolmogorov forward equation for the diffusion process {xp{t);t > 
0} on (0, 1) is 

86 1 d^ d 

(1.5) -^ = 2d^ ^^^^ - ^^'^> - 2^a^ ^^(^ - ^^"^^ ' 

with the initial condition (/)(p, x; 0) = 5{x — p). The solution is 

(1.6) <^(p, x; t) = 2(1 - r2)e'=(-^)e2- J^ ^^' ^"'"£' ^"'"^ pU^it), 

k=0 

where r = 1 — 2p, z = 1 — 2x, p^_,_2(i) := exp(— A^t), k = 0, 1,2, ... and — Ajt (0 < An < 



Kimura 



Ai < • • • ) are the eigenvalues of the generator. A plot of Aq is given in Figure 1 of 
(I1955I ). V-I^f^ {c,z) is the oblate spheroidal wave function (see Appendix): 

(1-7) vil\c,z) = Y,'fKc)Tnz), 

l>0 

3 
where T^{z) is the Gegenbauer function (may also be denoted by Cj^ (z)) and the summa- 
tion is over even values of / if A; is even, odd values of I if k is odd. Nik is the normalization 
constant of V-^^j^ {c,z). The probability mass at the exit boundaries are 

(1.8) f{p, 1; t) = 2(1 - r2)e^('^+i) f^ ^^^2^^^^|^f^^^^(l - P^.m 

k=o ^ ^^ 

and 

(1.9) /(p, 0; t) = 2(1 - r^)e^(-i) J] ^^ ^^^^ ^'' ^ 1 " Pl^^it)). 

In Section 2, the ancestral process {6„(t);t > 0} without absorbing states is studied. 
An explicit form of the probability distribution of hn{t) is obtained by using the moment 
duality between the ancestral process and the Wright-Fisher diffusion with directional 
selection. It is also observed that the distribution of the ancestral process converges to a 
stationary distribution. In Section 3, the rate of convergence is discussed. In contrast to 
the neutral process, the final rates of convergence are given by the largest eigenvalue for 
all the states. Bounds for the probability that 6„(t) is at the state 1 are obtained by an 



elementary martingale argument, which corresponds to the bounds (11.40 for the neutral 
process. In Section 4, the ancestral process with absorbing states is considered. It is shown 
that the first passage times of the ancestral process {bn{t);t > 0} at the states 1, 2, ..., n — 1 
are larger than that in the neutral process for all the states. By killing the modified process, 
in which the state 1 is the absorbing state, the form of the joint probability generating 
function of 6„(t) and the number of branching events is obtained. By using this formula, 
the expectation of the total length of the edges in the ancestral selection graph is obtained. 
In Section 5, the ancestral process of the whole population {boo{t);t > 0} is studied. It 
is shown that the process of fixation of the allele in the diffusion model corresponds to 
convergence of the ancestral process to its stationary measure. The time to fixation of an 
allele conditional on fixation is studied in terms of the ancestral process. It is shown that 
the density of the time to fixation of a single mutant gene conditional on fixation is given 
by the probability of the whole population being descended from a single real ancestral 
particle, regardless of the allelic type. In the neutral process, the density of the waiting 
time until the ancestral process hits the state 1 and the density of the conditional fixation 
time are given by the probability that the ancestral process is at the state 1. The property 
does not hold in the process with selection. 

2. Number of ancestral particles 

In this section, we will obtain an equation that relates the moments of the Wright- 
Fisher diffusion with directional selection to a Markov process that specifies the number 
of particles (real and virtual) that are present in the ancestral selection graph. To derive 
1 explo it the concept of duality from the theory of Markov processes 



th is result, we wi 



(lEthier and KurtzJ (|l98fil ): Section 4.4). If X = {Xt;t > 0} and Y = {Yt;t > 0} are 
Markov processes with values in sets Zx and Zy, respectively, then X and Y are said to 
be dual with respect to a function f{x,y) if the identity 

(2.1) KM{Xt,y)]=Ky[f{x,Yt)] 



7 
holds for every x G Zx and y € Zy- Duality is a useful concept because it allows us to 
use our knowledge of one process to learn about the other. Although there is no general 
procedure for identifying dual processes, duality can sometimes be deduced using simple 
generator calculations. Specifically, if Gx is the infinitesimal generator of the process X 
and Gy is the infinitesimal generator of the process Y , then the duality relationship shown 
in (12. ip will be satisfied if the identity 

(2.2) GJ{x,y) = Gyf{x,y) 

holds for all x and y. Here we think Gxf{x, y) as acting on the x-variable of the function 
f{x,y) for each fixed value of y. 

To apply these results to the Wright-Fisher diffusion with selection, it will be necessary 
to consider the frequency yq{t) = 1 — Xp{t){q = 1 — p) of the less fit allele, which is itself 
governed by a Wright-Fisher diffusion with generator: 

(2.3) Gyfiy) = ^y{l - y)f"{y) - 2cy(l - y)/'(y). 

Notice that the selection coefficient is negative in this case (c > 0). If we define the 
function f{y, n) = y", then a simple calculation shows that 

Gyfiy,n) = Q[/(y,n-l)-/(y,n)]+2cn[/(y,n+l)-/(y,n)] 

(2.4) = Gnf{y,n), 

where Gn is the operator defined by 

(2.5) Gnfin) = Q [fin - 1) - /(n)] + 2m[/(n + 1) - /(n)]. 

In other words, G„ is the infinitesimal generator of the birth and death process {6„(t); t > 
0} which keeps track of the number of ancestral par ticles in the ancestral sel ection graph. 



Ethier and Kurt J (jl98fil ) (Chapter 4, 



Because f{y,n) is bounded, we can use a result of 

Corollary 4.13) to deduce that the Wright-Fisher diffusion {yq{t);t > 0} and the birth and 

death process bn{t) are dual with respect to the function f{y,n): 



Theorem 2.1. 



(2.6) 



E[g^"W]=E[(y,(t)r], n = l,2,.. 



Because the right-hand side of this identity involves moments of the process {yq{t); t > 
0}, the process {bn{t);t > 0} is said to be a moment dual for {yg{t);t > 0}. The exis- 
tence of moment duals of Wright-Fisher diffusions with polynomial coefficients was first 
shown by IShigal (|l98ll ). and the explicit description of duality between the birth and 



death process and the Wr ight-Fisher diffusion with directional selection was discussed by 



Athreva and Swart 



2005l\. The 



graph (lAthreya and SwartI (J2005l )): 



duality can be proved in terms of the ancestral selection 



Proof. Partition an ancestral selection graph Q into disconnected subgraphs Qi,i = 1, 2, ... 
Let £t be the edges, or the ancestral particles, of a cross section of Q taken at time t 
backward. Then, bn{t) = \£t\. Each £q n Gi consists only of type A2 particles if and only 
if £t n Qi consists only of type A2 particles, since at least one type Ai particle survives 
from time t to if f^ Pi Gi contain type Ai particles. Here the ancestral selection graph is 
viewed forward in time. If a type Ai particle reaches a coalescence point, the number of 
type Ai particles increase by 1. If a type Ai particle reaches a branching point and meets 
another particle, the resulting particle is always type Ai. Thus, a sample consists only of 
type A2 particles if and only if the ancestral particles at time t consists only of type A2 
particles. D 

If a sample contains type Ai particles, then for n = 0, 1, ...; m = 1,2, ..., 
(2.7) E[(x,(t))™(y,(t))"] =E[(l-y,(t))™(y,(t))"] = j;(-l)^('^)E[g''-+"W], 

with the convention bo^t) = 0. In particular. 



j=0 



(2i 



E[xpit)] = E[l - y,it)] = E[l - g'-^W] = X;P[6i(t) = k]^, [ ■ Jp'q' 



k=l 



i=l 



9 

The expression follows immediately from the distribution of the number of particles in £f 
Since the ancestral selection graph is irreducible, a type Ai particle is sampled if and only 
if 8t contains at least one type Ai particle. The first summation is over \8t\ and the second 
summation is over the number of type Ai particles. 

By the Ito formula, we have a system of differential equations for the moments of yq (t) 



(2.9) -IJ^ = -{an+ Pn)^n + OinS,n-l+ Pn^n+l, n = l,2,... 



where ^„ = E[(yij(t))"]. The Kolmogorov backward equation for the ancestral process 
{bn{t);t > 0} without absorbing states is also given by ()2.9p . where ^„ = P[6„(t) = i]. 
Thus, (j2.9p with ^„ = E[g''"*-*''] holds. The isomorphism of these equations is a consequence 
of (1231). 

A realization of the ancestral selection graph consists of two disconnected subgraphs 
embedded in a diagram of a sample path of Xp{t) can be seen in Figure 1. The abscissa is 
the forward time interval (0, t) and ordinate is Xp{t). Thick lines represent the real geneal- 
ogy. Lines used by type A2 particles are dotted. The graph contributes to K[{xp{t))'^yq{t)] 
and 64 (t) = 4. If an ancestral selection graph consists of the upper subgraph only, it 
contributes to E[yg(t)] = E[g^i(*)] and bi{t) = 2. 



imi)) for / 



Using an integral transform by the Gegenbauer function ([Erdelvi et al. 
0,l,...;n = l,2,...;i = 0,l,..., 



(2.10) / Ti'{z){l + zni-zydz= '^^^''l-^^^^^^^^^^\ F2{-l,l + 3,i + l;2,i + n + 2;l), 
J-i (re + l)i+i 



where 3F2{—1, / + 3, i + 1; 2, i + n + 2; 1) is the generalized hyper geometric function, and 
with an identity (|5.2|) . it is possible to obtain an explicit expression for the probability 



10 

generating function of 6„(t), and we have 

E[g^"W] = E[y,itr]= [ {l-xr<Pix,p;t)dx + fip,0;t) 

Jo 



fe=0 -^^ I 



e'^'-l,... .2..cir-i,^n;.'ic,r)l^ ^^^ yll\^M pc^^^^t^ 



2A 
k=[) '" y. 

(2.11) = ^P[6„(t) = i]g\ 



i=l 



where 



{2cy ^ H),a + i),+2(i + i), 



r = 1 — 2p and /(p, 0; t) is the probabihty mass at given in (jl.9p . If n = oo, Ffcoo(c) = 
and /(p, 0;t) gives the probabihty generating function. Using a power series expansion in 
q of the Gegenbauer function 

(2-12) ^^(^)=(-i)T ^^'^;!|,Vir ^^' ^=°'''" 

we obtain an exphcit expression for the probability distribution of bn{t): 

(2.13) FMt) = 1] = vri + Se-^'^ f] ^^^^^^ |i^fcn(c) - ^^||^| Pfe+2(t), 
and 

(2.14) P[6„(t) = ^] = vr,+8e-2= f] ^ |f,„(c) - ^^^|^ | P^+2 (0 , ^ = 2, 3, ..., 



where 






2(i-l)! 



^ J-(j-l)! ,>~:^_i 2(z-j-l)!(i-j)! 
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and TTj are given in ()3.ip . Note that there are finite probabihties at the states n+1, n+2, .... 
The expected number of ancestral particles is 

(2.15) +8e-^^ E ^ E %^ l^Uc) - ^%^ I pU^it), 

i=2 k=0 ^f^ { ^ ) 

and the falling factorial moments are 

(2.16) 

nbnm] = ^!v^.e^^+8e-2^ f;b1. f] %^ |f,„(c) - ^^^j P^.+2(t), ^ = 2, 3, ... 

For small c, the probability distribution is approximately 
(2.17) 

n+l ^ r 

P[6„(t) = 1] = P[a„(t) = l]-2c+2c V(-1)'=(2A;-1) <^ | 

tl, l( 

and for i = 2, 3, ..., 






"+^ /' 1\fc-i 



[-lf-^{2k-l){i)k-l {k{k-l)[n]k- 



nhn{t)=^] = PK(t)=.] + 2c5,2-2cX:^^^^^^#^{ ,, 

ik"^ — k + 2i — 2)\n\h 1 o/n (^ — 1)j-iMj-i n , ■. ^, o^ 

(2-18) + ,/ . , ,,n , ■ o^ \ Pk{t) + 2c \. , ,/\^ P°-i(t) + O(c'), 



(fc - i + 1)(A; + i - 2)(n)fc j ^«^ ^ ' [i - l)!(n),_i 

with a convention [n]„+i = 0, where a„(t) is the number of ancestors of a sample of n 
neutral genes at the time t in the past. 

It is possible to obtain the solution of (|2.9p as a perturbation series in 2c, where the 
series is represented by eigenvalues of the neutral process. Let 

(2.19) in = d°^ + 20^(1) + (2c)2e(2) + . . . , n = 1, 2, ... 

Denote the rate matrix of the neutral process {a„(t);t > 0} by Qq = {qo,ij), where 
Qo,i+i,i = ctj+i^'Zo.M = —Oii for i = 1,2,... and other elements are zero. Let the Laplace 
transform of ^li (t) be f„ (A). It is straightforward to show that 

(2.20) i^(-'^ = {{Qo-XE)-'cyiXE-Qo)-'i'-'\o), i = l,2,... 
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where C = {cij) is given by cu = i,Ci^i+i = —i for i = 1,2,... and other elements are 
zero. Note that the inverse Laplace transform of the element in the n-th row and i-th 
column of the matrix {{Qq — XE)^^Cy{XE — Qo)~^ gives the j-th order coefficients in 2c 

0fF[bn{t)=i]. 

Let rn{t) be the number of branching events in the time interval (0,t) in the ancestral 
selection graph of a sample of n genes, where r„(0) = 0. The joint probability generating 
functions of 6„(t) and rn{t) satisfy a system of differential equation 



(2.21) 



dt 



-{an + V(in)in + On^n-l + v[3nin+l - (1 - t')/3nCri 



n = l,2. 



with the initial condition i^n(O) = q^, where ^„ = E[q^"(*''f'"" (*•*]. The solution is given by 
killing of the modified process {6„(t);t > 0} in which the selection coefficient is vc, and 
we have 



(2.22) 



^[q^r.{t)^r^{t)^ ^ ]g 



§''"(*) exp <{ -2c(l - v) I bn{u)du 



By setting r; = in the last expression, we obtain the identity 



(2.23) 



]E[<7''"^*\r„(t) =0] =E 



g""Wexpi-2c / an{u)du\ 



which shows the Poisson nature of the branching in the ancestral process {6„(t);t > 0}. 
The integral is the total length of the edges in the neutral genealogy without branching 
in (0,t). In particular, P[6i(t) = l,r(t) = 0] = e'^'^*. 



3. Convergence 



Standard results on birth and death processes (see, e.g., iKarlin and Tavloij ()l975l )) 
gives the stationary measure of the ancestral process {6n(t);t > 0}. It is straightforward 
to obtain the stationary measure 



(3.1) 



TTi := P[&„(oo) = i] 



(4c)^ 
i!(e4'= - 1) 



i = l,2 



•j ■"? '•'•> 
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which is the zero-truncated Poisson distribution. Since the ancestral process of the real 
particles is the neutral process {a„(t);t > 0} and the number of real particles becomes 1 
in finite time, vri is the probability that there are no virtual particles. 
It is clear from (f233|) and ^7\^ that for i = 1,2..., 

(3.2) vr, - F[bn{t) = i]= Oiplit)), t -^ oo. 

For small c, the final rates of convergence are approximately 

(3.3) 

hm {plit))-' {^1 - P[6„(t) = 1]} = 3e-24!i^ - ^" \ -2ce-^^5n,i + 0{c^), 

t^oo [n + l [n + l){n + 2)) 

(3.4) 

lim ipmr' {^2 - nbnit) = 2]} = -3e-"^{^ - Ml^zilEJ _ 2ce-2^5„,i + 0{c'), 
t^oo [n + 1 n + 1 J 

and 

(3.5) hm (p^(t))-i {vr, - P[6„(t) = i]} = _3,-2c ^ - j (^c)''^ ^ Q^^^i^i^^ • ^ 3^ 4^ __ 

In contrast to the neutral process, the final rates of convergence are given by the largest 
eigenvalue for all the states. In the neutral process, we have 

(3.6) hm (pO(t))-ip[a„(t) =i] = ^^, i = 1, 2, ..., n. 



The total variation norm has no simple form as in (jl.2p . 

A simple argument gives a bound for P[6„(t) = 1]. The event that the number of 
ancestral particles is 1 is a subset of the event that the number of real particle is 1, and 
we have 

(3.7) ¥[hn{t) = l]<W[an{t) = l], n=l,2,... 

An elementary argument on a martingale gives explicit bounds for F[6„(t) = 1], which 
corresponds to the bounds (jl.4p in the neutral process. Let r]{n] c) satisfy a recursion 

(3.8) (Ao -On -/?n)^(«;c) + a„r?(7i- l;c) +/?„r/(n + l;c) = 0, n = l,2,... 
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with the boundary condition r/(l;c) = —2c. Since rj is an eigenvector of the transition 
probabihty matrix of the ancestral pro cess ibn(t);t > 0), rjibnit ); c)( Po(t))~^ is a martin- 



gale to the ancestral process (see, e.g., 



Karlin and Tavloil (11975|)). Then, 



(3.9) 



K[r](bn{t);c)]=r]{n;c)pm- 



Although the explicit form of r/(n; c) is not available, it is possible to obtain an asymptotic 
form. Because of 



(3.10) 



Vin;c) 
r]{n — 1; c) 



1 + 



2An 



n{n — 1) 



+ 0(n-3), 



n -^ CX3, 



we deduce the asymptotic form r]{n; c) ?» r](oo; c)(l — 2Xo/n), where r/(oo; c) is a function 
of c. Although the explicit form of r/(oo;c) is not available, it is very close to 3exp(— c) 
(See Figure 2). 

For small c, 7/(n; c) can be expanded into a power series in c: 

(3.11) 

, , n-lf n'^ + n + 2 4(3n2 + lOn + 18) ^1 ^, o, 

T7 n; c) = 3 <^ 1 - -c + -^ -^c^ } + 0(c^ , n = 2, 3, ... 

'^ ' ^ n + ll (n-l)(n + 2) 25(n + 2)(n + 3) J 

Note that r/(oo;c) is not exactly equal to 3exp(— c). 
Lemma 3.1. ri{n;c) is monotonically increasing in n. 



Proof. Denote the rate matrix of the ancestral process {bn{t); t >0} hy Qc = iQc,ij), where 
Qc,i+i,i = ai+i)?c,M = — (oj + (3i),qc,i,i+i = Pi for i = 1,2, ... and other elements are zero. 
rj is an eigenvector of an oscillatory matrix E + Qc{2N)~^ which belongs to the second 
largest eigenvalue 1 — Ao(2A^)^^. An eigenvector of an oscillatory matrix which belongs 
to the second large st eigenvalue has exactly one variation of sign in the coordinates (see. 



Gantmached ()1959l ). pp. 105). Assume r]{i;c) > 0, i > L and rj{i;c) <0, 1 < i < L — 1. 
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Suppose I > L — 1. By an induction we deduce from (13. 8p that 

r]{l + l;c) ~r]{l;c) = —{??(/; c) - r/(/ - 1; c)} - Ao??(^; c) 

By taking t = oo in (|3.9p it follows that 

oo 

(3.12) J^ r?(i; c)^, = 0. 

i=l 

Thus, 

, oo 

(3.13) rj{l + 1; c) - r/(/; c) = -^^ V rj{v, c)T^i > 0. 

8^ ^'-1 ^tItl 



Next, suppose 2 < I < L — 2. We have 

(/-I) 



(3.14) ,y(; + i;c)-77(/;c) 



^ I v{2; c) - r/(l; c) - ^ ^ ??(i; c)7r, I > 0. 



(4c)' 
Finally, r/(2; c) - r/(l; c) = Aq > 0. D 

From Lemma 13.11 it follows that 

P[6„(t) = 1]7?(1; c) + P[6„(t) > l]r?(2; c) < E[r?(6„(t); c)] 
(3.15) < P[6„(t) = l]77(l;c)+P[6„(t)>l]77(oo;c), 

here we note that there are finite probabilities at the states larger than n. Then, from 
(|3.9p we have the following bounds: 



Theorem 3.2. // rj{n; c) satisfies the recursion 113. 8\) . then 

r?(n; c)e"^o* + 2c „ ,, , r?(n: c)e"^o* + 2c 

3.16 '^ ' <i-^[h^(t) = l]<J^^±- , n = l,2,... 

??(oo; c) + 2c Ao 

Figure 3 shows the bounds when c = 0, 0.1, 1 and 8 for n = 10. When c > 0, in contrast 
to the neutral case, the states larger than 1 have finite probabilities in the stationary 
measure (13. ip and the upper and lower bounds do not converge to the same value. Figure 
3b shows that the upper bound is sharp for small values of c and loose for intermediate 
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values of c (Figure 3c), and that both the upper and lower bounds converge to 1 as c tends 
to infinity. 



Let 



4. First passage times 



(4.1) 



Wli:=mi{t>0-bnit)=i}, 



1,2, 



and {bli{t);W^i > t > 0} be a modified process, where there is an absorbing state at 
1, or the ultimate ancestor. The modified process is the same as that introduced for 



ancestral recombination graph, where 4c is rep 



(|Griffithsl dlHail)). Theorems 1, 2 3 in ^Griffiths 



modified process was studied by 



aced by the recombination parameter p 



199lh hold fo r the modified process. The 



Krone and Neuhauseij (jl997l ). Here, modified processes 



{6^(t); ly^ j > t > 0}, where there is an absorbing state at i = 1,2, ...,n — 1, are studied 
to discuss the first passage times of the ancestral process {6„(t);t > 0} at the states 



l,2,...,n-l. 



It is possible to show that the expected first passage times of the ancestral process 
{bn{t);t > 0} at the states 1,2, ...,n — 1 are larger than th ose in the neutral process 



{an{t);t > 0}. E[VF^i] is given in 



Krone and Neuhauseii (|l997l ) 



Theorem 4.1. Let 



(4.2) 



Wli := inf{t > 0; a„(i) = i}, i = 1, 2, 



Then, 



(4.3) 



n— 1 oo 



nwi,] = 2^Yl 



i4cy 



where W^ i is the time to the ultimate ancestor. 



>nwli], i = l,2,...,n-l. 
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Proof. The theorem follo ws from standard results on birth and death processes (see, e.g., 



Karlin and Tavloij (|l975l )). The modified processes {bnit);W^^ > t > 0} hit the states 



i = 1,2, ..., n — 1 in finite time with probability one, since 



(4.4) 



V TT ^ = ^ ' V 



ml 



Pk {i - 1)! f^. (4c)^ 



oo, 



i = 1,2, ...,n — 1. 



m=i k=i+l 

From the Kolmogorov backward equation for the modified process {Vn{t);W^^ > i > 0}, 
which is (12. 9p for n = i + l,i + 2,... with ^„ = P[6^(t) = i] and the boundary condition 
S,i = 5{t), the expected first passage times satisfy a recursion for i = 1, 2, ..., n — 1 



(4.5) {an + f3n)C{n)-anC{n-l)-PnC{n + l) = l, n = i + l,i + 2, ...,n -2 

with the boundary condition C,{i) = 0, where C,{n) = E[W^ J. It is straightforward to solve 
the recursion and obtain 

oo n—2 j+1 oo n— 1 oo 



(4.6) nwi.] = Y.im + Y. n i^ E 7. = 2EE 






k=i j=0 



(4c)^ ^ 

(^)i+2 ' 



i = l,2,...,n-2. 



and 

(4.7) 
where 

7i 



E[TV,Vi] = E^- = 2E 



j=0 



1 



Oj+i i(i + 1) 



2 /?i+i/3i+2 • • • A 
, 7m ~ 



(4c)^' 

(^)j+2 ' 



2(4c) 



aj+l«i+2 • • • "mam+l {i)m~i+2 



m = i + l,i + 2, ... 



It is clear from (14. 6p and (I4.7P that 



(4.8) 



nWli>2"t^^=2(^\-l.)=ElWl]. 



z = 1, 2, ...,n — 1. 



L, .^, ..., 



D 



As c — > oo, for i = 1,2, ..., n — 1, 



n-l 



(4.9) E[I^„y=2^ 



fc+2 ) 



k=i 



^.k{k + 






k+2 



n-1 -^£- 

Eefc+2 
, .k{k + l) 



OO. 
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Corollary 4.2. For the whole population (n = ooj, the expected first passage times are 



(4.10) ^[^^.] = 2E(-Ti4;7' ^ = i'2'- 

Proof. It follows immediately from an identity 



°° 1 1 

fr^(fc)j+2 (j + l)Wj+i 



D 



It is straightforward to obtain higher moments of the first passage times of the ancestral 
process {hn{t)]t > 0} at the states 1, 2, ..., n — 1 in the same manner. The second moments 
E[(W^ j)^] satisfy a recursion 

(4.12) {an + /3n)C(«) - a„C(n - 1) - /3„C(n + 1) = 2E[I^,^J, n = i + 1, i + 2, ..., n - 2 

with the boundary condition C,{i) = 0, where C(^) = ^I^[(^ni)^]- However, there is no 
simple form for the density as in p.3p . The Laplace transform of the first passage times 
of the ancestral process satisfy a recursion for z = l,2,...,n — 1 

(4.13) (A + a„ + /3„)C(n) - QnC(n - 1) - /3nC(n + 1) = 0, n = i + l,i + 2, ...,n - 2 

with the boundary condition C,{i) = 1, where C(^) = lE[e^ "'']. 

The joint probability generating function of h\{t) and r„(t) satisfies a system of differ- 
ential equation (|2.21|) with ^„ = E[g^"(*)f *""(*)]. By taking t = oo, we have 

(4.14) Q = -{an + Pn)in+anin-l + vPnin+l^ n = l,2,.., 

with the boundary condition ^i = 1, where ^n = E[i;''"'°^']. The form of the probability 
generating function of r(oo) is 



(4.15) E[i;''"(°°)] = E 



exp < — 2c(l — v) I h^{u)du > 
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while the explicit fo r m of the probability generating function is given by Theorem 5.1 in 



Ethier and Griffitha (|l99d ). where p is replaced by 4c, and we have 

RJv) 



(4.16) 
where 



Ergr„(oo)i 



Ri{vy 



Rn{v) 



Jo 

(n-1)! 



(4c(l-«))„ 

(n + i - l)!(-4cf) 



Fi (n; 4c(l — v) + n; —Acv) 



. „ (4c(l 



i=0 



{4c{l-v)+n)n+iil' 



(|4.15|) provides a way to compute the expectation of the total length of the edges in the 
ancestral selection graph in the time interval (0, W^ i), and we have 



(4.17) 



E 



ws 



hn{u)du 



-. oo n— 1 -. 

-E[r„(oo)] = Y.(^cf~' E J- 



fe=i 



m=l 



m k 



It is possible to obtain the probability that the modified process {6^(t); W^^ > t > 1} 
hits the states n + 1, n + 2, ... 



j = 3,4, 



Theorem 4.3. Let z{l) = 0, z{2) = I, and 

(4.18) .(j)^l + f + ^ + ... + °7'-7-' =gji^. 

P2 a2P3 P2P3 • • • Pj-1 ^ (4c)'' 

Then, the probability that the modified process {bl^{t);W^i > t > 1} hits the states m 
n + l,n + 2, ... is 

z{n) 



(4.19) 

Proof. The theore r a fol 



lP[W^n,l > W'n 



z{m) 



ows from standard results on birth and death processes (see, 
Karlin and Tavloij ( 19751 ). pp. 323). It is straightforward to show that z{h\{t)) is a mar- 
tingale for the modified process. minjW^^, W^^} is a Markov time with respect to the 
modified process. We apply the optional sampling theorem to conclude that 



(4.20) z{n) = E[6i (min{W„^,i, ly^^})] = ^[^n,i > W^^Jz{m), m = n + 1, n + 2, 



20 

D 



Remark 4.4. For small c, ¥[Wf^i < W,^^] can be expanded into a power series in c. 



(4.21) 



Wn,l > W'n, 



(4c) 



[m — 2\r 



+ 0{c 



,m— n+l> 



m = n + l,n + 2, , 



Figure 4 shows the hitting probabiHty for n = 10 and 50 as a function of c. It can be 
seen that the hitting probabihty grows quickly around a critical value of c. For n = 50 
and m = 51, the probability grows linearly until c is smaller than 4.5 (See Remark 4.4). 
Then it approaches 1 quickly. 



5. TIME TO FIXATION 

In studying evolutionary processes from the standpoint of population genetics, the prob- 
ability and the time to fixation of a mutant gene play important roles. The expected time 



to fix ation of a muta nt gene condi tional on fixation was obtained by 



(119691 ). Furthermore, 



Ewena (jl973l ) and 



Kimura and Qhta 



Maruyama and Kimura (1974 ) showed that the 



expected length of time which it takes for an allele to increase frequency from qio y {> q) 
on the way to fixation is equal to the expected length of time which the same allele takes 
when its frequency decrease from y to g on the way to extinction. The time-reversibility 

soiourn time does 



Maruvamal (|l972l ). 



property is equivalent to the property that the density of the expectec. 

not depend on the sign of the selection coefficient, which was shown by 

While these results are well known, their interpretation in terms of the ancestral process 

of the whole population {boo{t)',t > 0} are interesting. 



The fixation probability was obtained by solving the Kolmogorov backward 



e quati on 



for the Wright-Fisher diffusion with directional selection {xp{t);t > 0} (JKimura 
The fixation probability of the allele Ai is 



mm) 



(5.1) 



ui{p) 



-icp 



-Ac 
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and the fixation probability of tlie allele ^2 is 1 — ui{p)- It follows from ()1.9p that 



00 tAi) I 



.(1), 



(5.2) 



1 - uiip) = 2(1 - r2)e'=(^-^) J2 ^'' ^"'''^^'' ^""'^ 



fe=o 



2Xk 



It is possible to obtain the fixation probability from the stationary measure of the ancestral 
process (|3.ip . If the allele A2 fixes in a population, the ancestral particles of the whole 
population in infinite time backwards consist of type A2 particles only, and we have 



(5.3) 



E[q' 



,600(00)] 



X^vr^g^ 



i=l 



=4cg _ -^ 



1 -ui{p). 



The relationship between the fixation probability and the pro bability generating functio n 



Athreva and SwartI (|2005l ). 



of the stationary measure is a special case of Theorem 2 (f ) in 

The density of time to fixation of the allele A2 conditional on fixation has a genealogical 
interpretation. Let 



(5.4) 



ro^:=inf{t>0;y,(t) = l}. 



Then, it follows from the expression 



(5.5) 



p[ro" < t\T^ < 00] 






that the process of fixation of the allele in a diffusion model, in which the left hand side 
converges to one as t — > cxd, corresponds to convergence of the distribution of the ancestral 
process ¥[boo{t) = i] to its stationary measure ttj as t ^ 00. 



The expected time to fixation of the alle^ 
solving the Kolmogorov backward equation (jKimura and Ohta 
and 



e Ay, conditional on fixa ti on was obtained 



ilMSl), 



Maruvamal(ll972l )). 



by 



(5.6) 



EiTo^lTo^ < 00] = / ^q,y)dy, 
Jo 
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where ^{q,y) is the density of the expected sojourn time of the ahele A2 at frequency y 
in the path starting from frequency q and going to fixation, and 



^iq,y) 



y>q, 



(5.7) 



S{y)S{l-y) 
2c2/(l-y)5(l)' 

S{y) iS{l-y) S{q - y) 

S{q) 



2cy{l-y) \ S{1) 
and S{y) = exp(4cy) — 1. Then, 

1 S{y)S{l-y) 



y <q, 



(5.8) 
where 



E[TS\TS < 00] 



dy 



Siy)Sil-y) 
2cy{l-y)S{l] 



dy 



2cy(l - y)5(i: 



TTl ^;^^;^ (4c 



^ Siy)Siq-y) 
2cyil-y)Siq) 



dy, 



e^2 2^Z^ 



\i+j /-l 



8c2Z_.^ i!^-, ^^ 



f-\l-yy-'dy 



2cA.(fc + i)!Z.i(A:_, + i) 
,^^gi^..r(4c)^ 



fc=0 



(fc + 2)! ' 



Hk = l + l/2 + --- + 1/k, and 
'•^ 5(y)5(g-y) 



2cyil-y)Siq) 



dy 



+j riyi-\q-yy 






^cS{q)f^^P^ i\3\ Jo 

1 °2, °2, (4cg)»+i 
2c5((?)^^^i(z + i)! 



y 



dy 



2Fi{l,i,i + j + l;q) 



1 ^^(4cg)^+^-^ {^,q' 






It is possible to obtain the expected time to fixation of the allele A2 conditional on fixation 
from the distribution 600 (0; ^^^ ^^ have 

1 



E[TS\TS < 00] 



,00 ^ 
Ei^l '^iQ' ./n dt 





00 



i=l 



(5.9) 



E- 



=i^«9 ~^ Jo 



Y^{F[b^{t) = i] - 7ri}q' 
(t) =i]- 7ri}dt. 



dt 
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From the two expressions (15. 8p and ()5.9p . an identity at g = follows immediately. 

Hk+ii'ic)'' 



(5.10) 






'"'-IT. 



{k + 2)\ 



k=0 f^ fc=0 

It is straightforward to obtain similar identities by comparing (|5.8p and (|5.9p in each 
power of g. Moreover, explicit expression for the higher mo ments of the time to fixa- 



19721 ) ■ and they produce similar 



tion, conditional on fixation, are available ([Maruvamal . 
identities. 

The density of the time to fixation of a single mutant gene conditional on fixation has 
interesting properties. Let 



(5.11) 



Tf :=inf{t >0;xp(t) = 1}. 



" ?hen, from a time-reversibility argu ment on the conditional diffusion process (JEwens 



m7± 



Maruyama and Kimural (Il974l )). we have 



(5.12) 



lim ¥m < t\n < oo] = lim P[rf < t|rf < oo] 



At) = 1] 



VTi 



The same density holds for a mutant gene of allele Ai and a mutant gene of allele A2- 
This property has an intuitive genealogical interpretation. The conditional density is given 
by the probability of the whole population being descended from a single real ancestral 
particle. Since there is no variation in the population, selection cannot have an effect on 
it and consequently, the conditional density should not depend on the allelic type. 

Theorem 13.21 gives bounds for the density of the time to fixation of a single mutant gene 
conditional on the fixation, and 

(5.13) r < lim FlTn < tlTn < 00 < r — — . 

^ vTi TTiAo q-*o ^ ^ ' ^ vTi 7ri(?7(oo;c) + 2c) 

The bounds are useless when c is large, since the upper and lower bounds do not converge 
as c — > 00. When c is small, the lower bound is sharp. Figure 5a shows the bounds 
when c = 0.1. Figure 5b gives the conditional density and the derivative of the lower 
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bound 7ri^^r/(oo; c) exp(— Aot). It can be seen that the tail of the conditional density is 
well characterised by the largest eigenvalue. 

It is worth noting that the identity ()5.12p gives the following identity in the distribution 
bn{t). Its interpretation in terms of the ancestral process {bn{t); t > 0} is unclear. 

Remark 5.1. i5.13\) gives 

(5.14) P[6oo(t) = 1] = lim e-''^y2(-^)'^'(l)^Mt)]- 

Proof. (|5.12p is equivalent to 

(5.15) ,_^/(p^^,__^/(p,Oa)_P[Mt) = ll 



+0 U\{p) 9^0 1 — U\{j)) VTi 



where 



lim — - — = hm iim — = lim lini --\ 



nq'"'^''>] .^ „-4c^, ,,k-,ifn\nbkit)] 



(5.16) = lim lim ^ (-1)'= f'^'j ^li^lH = lim e'^-^V (-1)'=+^ f 



kj TTi 

D 



In the neutral Wright-Fisher diffusion, the density of time to fixation of a mutant gene 
conditional on fixation follows 

(5.17) limP[ro° < t|ro° < oo] = P[aoo(t) = 1], 

where Tq is the time to fixation of a mutant gene in the neutral Wright-Fisher diffusion. 
From (j5.8p . the expected time to fixation of a mutant gene conditional on fixation has a 
simple form 



(5.18) lim EfTo^lTo^ < oo] = 4^i f] ^^^ < l^mM < ^] = 2, 

where the inequality holds from the following lemma: 
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Lemma 5.2. The density of expected sojourn time of the allele A2 at frequency y in the 
path starting from frequency and going to fixation satisfies 

^'■''^ 2cyil-y)Sil) < '' '<'<'■ 

Proof. The inequality is equivalent to 

p4cj/ _ 1 p4c(l-j/) _ 1 

(5.20) < Ac{e^^ - 1), 

y i-y 

or 

r,9n f.^ (4c)V(i-^r^- ^ {Acf 

^ ' ^^{j + iy.{i-j + iy. ^{i + i)V 

The inequality follows from an inequality 

^Q(j + l)!(i-j + l)! (^ + l)!^t^ J!(^-J)! (^ + 1)!' 

D 

As c becomes large, P[6oo(0 = 1] decreases, while the expected fixation time of a mutant 
gene conditional on fixation decreases. It is straightforward to show that the inequality 
for the expected fixation time (jS.lSp is equivalent to an inequality 

(5.23) /'"|^M!l^_P[a^(i) = l]|dt>0. 

In the neutral process, both the density of the waiting time until the ancestral process hits 
the state 1 and the density of the conditional fixation time are given by the probability 
that the ancestral process is at the state 1 (jl.3l5.17p . It follows that 

/>oo 

(5.24) E[W<^ J = lim E[ro°|ro° < 00] = / {P[aoo(t) = 1] - \\dt = 2. 

9^0 Jo 



In contrast, in the processes with directional selection, we have 

(4c)^' 
(j + l)(i + l)! 



(^•25) E[^^J = 2i;7--W7TW>2' 



j=0 



26 

while 

(5.20) linjEKIT.^ < c„, ^ /°° {'-MH^ - l| ., ^ 4.. f ^^^ < 2. 

•J^o Jo I ^1 J ^ U + 2)! 

6. DISCUSSION 

In this article, the ancestral process {6„(t); t > 0}, specifying the number of branches in 
the ancestral selection graph, was investigated by exploiting the moment duality between 
this process and the Wright-Fisher diffusion with directional selection. An explicit formula 
for the probability distribution of 6„(t) was derived. Although this expression cannot be 
given in closed form, since it involves eigenvalues and coefficients which are determined 
by an intractable three-term recursion relation, it is possible to expand the probability 
distribution as a perturbation series in 2c. This expression is given in a closed form for 
each order of the perturbation and is accurate when c is small. Bounds for the probability 
that 6„(t) is at the state 1 is obtained. When c is small, one of the bounds is sharp. The 
density of time to fixation of a single mutant gene conditional on fixation was shown to 
be given by the probability of the whole population being descended from a single real 
particle, regardless of the allelic type. Thus, the bounds for the probability that bn{t) is 
at the state 1 give bounds for the density of the conditional hitting time. It was shown 
that the tail of the conditional hitting time is well characterised by the largest eigenvalue. 
The probability that the process hits states larger than the initial state before the process 
hits the state 1 was obtained. According to the formula, the number of branches in the 
ancestral selection graph grows rapidly when c is larger than a critical value. One of the 
difficulties of simulating the ancestral selection graph is keeping track of large numbers of 
branches when c is large. Specifically, when c is larger than the critical value, enormous 
number of branches emerge and it will be difficult to simulate the ancestral selection graph. 

If a sample consists only of type A2 particles, the probability distribution of the ancestral 
particles, all of which are A2, is bn{t) (see Theorem 2.1). If a sample contains type 
Ai particles, the joint probability distribution of the number of the Ai particles and 
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the number of the A2 particles is interesting. However, it seems that the expression 
of the moments in the Wright-Fisher diffusion (j2.7p does not give any insights into the 
joint probabihty distribution, except in the case when a sample consists of a single Ai 
particle. The time-reversibility argument of the conditional diffusion process gives an 
identity, whose interpretation in terms of the ancestral process is unclear (see Remark 
5.1). The interpretation of the time-reversibility in terms of the ancestral process needs 
further investigation. 

In this article, the fixation process of a mutant gene in the Wright-Fisher diffusion with 
directional selection, which has been well studied, was interpreted in terms of the conver- 
gence of the ancestral process to its stationary measure. This approach might be useful 
for other population genetical models for which less is known about the fixation process. 
Examples include diffusion processes with frequency-d ependent selection, with m ultiple 



selected alleles, and of spatially-structured populations (JAthreya and Sward . 



20051). 



Appendix A. The oblate spheroidal wave function 



Th e oblate sphero i dal wa ve function ^^^ (c, z) can be represented by expansions of the 



form (Stratton et al 



(1941 



)) 



(A.l) 



vll\c,z)=Y,'ftic)Tnz), k = 0,l,... 



l>0 



'his notation was used in 



Kimura 



(Il94lh and {1 - z'^p Sik+i{c, z) in 



195511. It was d enoted by V^f, {—ic, z) in 



(1)^ 



Stratton et al 



Flammerj (119571 ). From the orthogonal properties of the 



Gegenbauer function it is shown that 



(A.2) 



.2at/(1) 



(1)/ 



(1 - z')Vl^>{c, z)V^i'>{c, z)dz = 5k,iNik 



where 



'^^•^^^E'^^M^iftie)?- 



l>0 



(2/ + 3) 
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Note that 



(A.3) V}l\c,l) = i^'(/ + !)(/ + 2)/f(c), V}l\c,-1) = (-l)'=F«(c,l). 

l>0 

The coefficients fj^ic) satisfy a three-term recursion in the form 
(A.4) A+2/f+2(c) + Biftic) + Ci,2fl2{c) = 0, 

where 

(/ + !)(/ + 2) ^ l{l + 'i)-bk 2P + QI + 1 ^ {l + l){l + 2) 



(2/ + l)(2/ + 3)' ' c2 (2/ + l)(2/ + 5)' ' (2? + 3)(2/ + 5)^ 

and 6jt = 2Ajt — 2 — (?. fi{c) = for odd Hf fc is even and for even / if fc is odd. ()A.4p can 
be developed as a continued fraction. 

/; M+2 C1-2A1 C2A4 A2 



ft+2 Bl- ^'-2- B2- Bo 



1 = 0,2,... 



(A.5) _ Ai^2Q.2Ai qsA^A, ^^i^3^__ 

and 

(A.6) % = C^ ^i+^Ci+2 __ / = o,l,.. 

^ ^ ft Bi+2- Bi+,- 

bk is determined by the condition that the reciprocal of the ratio fi/ fi+2 by (jA.SP must 
equal the value of fi+2/ fi obtained from ()A.6p . Then, the continued fractions provide a 
way to compute arbitrary coefficient. 

For small c, the eigenvalue can be expanded into a power series in c. 

.,^. . (fc + l)(fc + 2) , (fc + l)(fc + 2) 2 , ^, 4^ 

If we set /fc (c) = 1, then 

(A.8) 
fk ( -. _ {k + l){k + 2) o 4 fc |-„^ _ (A: + l)(A: + 2) ^ , -n/ 4^ 



and other coefficients are zero up to 0{c' 



A\ 
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Figure 1. A realization of the ancestral selection graph embedded in a diagram of a sample 

path of Xp{t). 

Figure 2. 77(00; c) (dots) and 3exp(— c) (line). 

Figure 3. 1 — P[6io(i) = 1] (line) and the upper and the lower bounds given by Theorem 

32] (dotted lines). 

Figure 4. The probability that the ancestral process {bnit);t > 0} hits states m{> n) 

before the process hits state 1 (ultimate ancestor). 

Figure 5. (a) The cumulative probability of the density of time to fixation of a single 

mutant gene conditional on the fixation (j5.12p (line) and the upper and the lower bounds 

given by (jS.lSp (dotted lines), (b) The density of the conditional fixation time (line) and 

vrj" 77(00; c) exp(—Aot) (dotted line). 
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Figure 2. 
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